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1 Introduction 



A pattern formed by a continuum of periodic orbits is the usual look of phase space for 
a conservative planar system. Take, for instance, the harmonic oscillator x + x = 0. The 
set of circles x 2 + x 2 = a 2 , with a a positive real number representing the amplitude of 
the oscillation, fills the whole plane (x,x). Nevertheless, dissipation and amplification are 
always present in real dynamics, since the oscillating system is usually embedded in an 
interacting external medium. The effect of this external medium is implemented in the 
evolution equations of the oscillator, with more or less accuracy, through a nonlinear term 
h(x, x) that depends on the position and velocity of the oscillator. In general, even if 
this new term is only a slight perturbation governed by the small real parameter e, the 
coexistence of the infinitely many periodic motions of the original unperturbed system 
is destroyed and only a finite number of them survives in the new context given by the 
equation x + x + eh(x, x) = 0. Dissipation and amplification are balanced on those orbits. 
These localized periodic motions that verify strict energetic balance conditions are called 
self-sustained oscillations or limit cycles 0121 • 

Lienard equations, 



with h(x,x) = f(x)x and f(x) any real function, model the dynamics of specific planar 
systems where limit cycles can be found. When f(x) is a polynomial of degree N = 2n + l 
or 2n, with n a natural number, Lins, Melo and Pugh have conjectured (LMP-conjecture) 
that the maximum number of limit cycles allowed is just n ||. It is true if iV = 2, or iV = 3 
or if f(x) is even and N = 4 [31 E]. Also, it has been claimed its truth in the strongly 
nonlinear regime (e — * oo) when f(x) is an even polynomial [3]. There are no general 
results about the limit cycles when f(x) is a polynomial of degree greater than 5 neither, 
in general, when f(x) is an arbitrary real function [5]. When fix) = x 2 — 1 we have the 
van der Pol oscillator, x + e(x 2 — l)x + x = 0, wich is a particular example of Lienard 
system that has been very well studied. It displays a limit cycle whose uniqueness and 
non-algebraicity has been shown for the whole range of the parameter e [Jj- Its behavior 
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runs from near-harmonic oscillations when e — > to relaxation oscillations when e — > oo, 
making it a good model for many practical situations 8 . 

Different non-perturbative approaches that allow to obtain information about the num- 
ber of limit cycles and their location in phase space have been proposed in the last years 
when f(x) is an even polynomial. A method that gives a sequence of algebraic approxi- 
mations to the equation of each limit cycle can be found in [§], and a variational method 
showing that limit cycles correspond to relative extrema of certain functionals is explained 
in 0. 

In this paper, we are interested in determining the number, amplitude and shape of 
the limit cycles emerging in the weakly nonlinear regime (e — > 0) of Lienard equations 
by standard perturbation techniques. This work completes a previous work jS] in which 
these questions were successfully answered and solved in the strongly nonlinear regime 
(e — ► oo). In fact, the algorithm here explained is an alternative to a more general, 
although technically more sophisticated, perturbative method proposed by Giacomini et 
al. jlUl II 1 j on the same question. Thus, we propose a new algorithm which permits an 
easy computation up to an arbitrary order 0{e ) of the amplitude a and form y(x) of the 
limit cycles. Moreover, this simple algorithm permits to show some interesting properties 
of symmetry of the approximate solutions. The article is organized as follows. A few 
symmetries of the Lienard equation are recalled in Section 2. The method is detailed 
in Section 3. Some illustrative examples are given in Section 4. Finally, we present our 
conclusions. 



2 Symmetries and an integral equation 

In order to study the limit cycles of equation Q with the time variable being implicit, it 
is convenient to rewrite it in the coordinates (x, x) = (x, y) in the plane, with x(t) = y(x) 
and x(t) = y(x)y'(x) (where y'(x) = dy/dx): 

yy' + ef{x)y + x = 0. (2) 
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A limit cycle C\ = (x,y±(x)) of equation @ has a positive branch y+(x) > and a 
negative branch y~(x) < 0. They cut the x-axis in two points (— a_,0) and (a+,0) with 
a_,a+ > because the origin (0,0) is the only fixed point of Eq. @. Then every limit 
cycle C\ solution of Eq. ((2j) encloses the origin and the oscillation x runs in the interval 
— a_ < x < a+. The amplitudes of oscillation a_,a + identify the limit cycle. The result is 
a nested set of closed curves that defines the qualitative distribution of the integral curves 
in the plane (x,y). The stability of the limit cycles is alternate. For a given stable limit 
cycle, the two neighboring limit cycles, the closest one in its interior and the closest one 
in its exterior, are unstable, and vice versa. 

When f(x) is an even function, the inversion symmetry (x,y) <-> (— x, —y) of Eq. (j2j) 
implies y+(x) = — y~(— x) and then a\ = 0,2 = a. Therefore we can restrict ourselves to the 
positive branches of the limit cycles with —a < x < a. In this case, the ampli- 

tude a identifies the limit cycle. The parameter inversion symmetry (e, x, y) «-> (— e, x, —y) 
implies that if C\ = (x,y±(x)) is a limit cycle for a given e, then C\ = (x, —y^(x)) is a 
limit cycle for — e. In consequence, the amplitude a of the limit cycles in Lienard systems 
is an even function of e. Moreover if C\ is stable (or unstable) then Ci is unstable (or 
stable, respectively). Therefore it is enough to consider the limit cycles when e > for 
obtaining all of the periodic solutions. (The limit cycles for a given e < are obtained 
from a reflection over the x-axis of those limit cycles obtained for e > 0). 

Another property of a limit cycle can be derived from the fact that the mechanical 
energy E = (x 2 + y 2 )/2 is conserved in an oscillation: 



Taking into account that dE = {yy' + x) dx = —ef(x)ydx, if equation © is integrated 
along the positive branch of a limit cycle, between the maximal amplitudes of oscillation, 
we obtain: 



The solutions y+{x) of Eqs. and |2J) and vanishing in the extremes, constitute the 
finite set of limit cycles of Eq. Q. 





(3) 
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3 The perturbative method 



For later convenience, we rewrite the Lienard equation (J2J with a different notation: 

uv! + ef(t)u + 1 = 0, (4) 

where f{t) is an even function of the real variable t and u' means derivative with respect 
to t. From the discussion of the above section we know that the amplitudes of the limit 
cycles are symmetric: if u{t) is a limit cycle of amplitude a then u(—a) = u(a) = 0. After 
the following change of dependent and independent variables: 

t = ax, u = ay, a > 0, (5) 

equation (J1J) reads 

yy + tf(ax)y + x = 0, (6) 

where a is now an explicit positive parameter of the equation, and y' means now derivative 
with respect to x. Every limit cycle u(t) of (JU of amplitude a is transformed into a limit 
cycle solution y(x) of © of amplitude 1 verifying y(— 1) = y(l) = 0. The main result 
of this paper is to establish in these new variables (x, y) a recursive algorithm capable of 
approximating a limit cycle solution y(x) and its amplitude a in a power series of e up 
to an arbitrary order. The approximated limit cycles in the original variables (u, t) are 
obtained by undoing the change of variables © at each order 0(e N ) of approximation. 
In order to show our results we need some previous definitions. 

Definition 1. For N = 0, 1, 2, we denote by y( N \x) the function that approximates a 
limit cycle solution y(x) of Q up to the order 0(e N ): 

N 

^)(x)e^ € \(x). (7) 

ra=0 

We denote by the approximation to its amplitude a up to the order (^(e^): 

N-l 

71=0 

The coefficient functions y n {x) and the coefficients a n are computed by means of the 
algorithm given below. 
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Definition 2. For n = 0, 1, 2, TV, we denote by f n (ao,ai, ...,a n ;x) the coefficients of 

the formal Taylor expansion of the function g(e) = f(a^x) at e = 0: 

/ n \ N 
flxJ2 ««6 n = £ / n (oo, ai, a n ; x)e n + ©(e^ 1 ). 

V n=0 / n=0 

The first few coefficients / n (ao,ai, ...,a n ;x) are: 

/o(a ;x) = f{a x), 

fi(ao,ai;x) =x/'(a x)oi, 
/ 2 (a ,ai,a 2 ;x) = xf'(a x)a 2 + -x 2 f" (a G x)a\, 
f 3 (a Q ,ai,a2,a 3 ;x) = xf'(a x)a 3 + x 2 f" (a x)a 1 a 2 + ^/'"(ao^a?. 

The general form of the coefficients f n (ao, a±, a n ; x) is given in the following lemma 
whose proof is straightforward and we do not reproduce here: 

Lemma 1. For n = we have /o(ao! x ) = f(a>ox) anc b for n = 1, 2, 

/ n (a ,ai,...,a n ;x) = £ 6„, m (ai, a n )x m / (m) (a x), (9) 

m=l 

where 

<t(1) cr(2) cr(n) 

In these formulas n > m > 1 and the above sum is performed over the following subset 
S n ,m of the set of all of the possible mappings a between the sets {1, 2, 3, n} and 
{0,l,2,3,-,m}: 

{n n \ 

a : {1,2,3, ...,n} -> {0,1,2,3, ...,m}; £ cr(/c) = m, ^#) = n . 
fc=i fc=i J 

We also define 6o,o = 1 and 6 nj o = for n > 0. 

Definition 3. Suppose that yo(x), yi(x),..., y n (x) are given. For n = 0, 1,2, we define 
the functions: 

n r l 



/3 n (a ,ai, ...,a n ) = V / / m (a , ai, a m ; x)y n - m (x)dx = 

£ £W / x k y n - m (x)fW(a x)dx. (10) 



m=0 fc=0 
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Algorithm. We compute the coefficients (y n , a n ), n = 0, 1, 2, N , of expansions Q) and 
in the following way. 



Step 1; Setyo(x) = \fl — x 2 and take one of the solutions ao of the (in general nonlinear) 
equation /?o(ao) = 0: 

A)(°o) = j fo(a x)y (x)dx = 0. (11) 



i 



Step 2: For each n, with n = 1,2, N , compute, firstly, y n {x) by means of the formula 



Vn{x) 



( j n—1 n—1 



2/o 0) 



- X y m (a;)yn-m(a;) + X / yn-m-i(t)f m (ao,ai, -,a m ;t)dt> = 

m=l m=0 ^ J 

re— 1 n—1 m „ x ~\ 

X y m (a;)yn-m(a;) + X X W / t k y n ^ m -i(t)f^(a t)dt \ (12) 

71=1 m=Ok=0 1 J 



yote) [2 

and, secondly, solve for a n the linear equation 



f3 n (a Q ,ai,...,a n ) = 0. (13) 



The derivation of this algorithm is given in the following theorem. Its detailed form 
for N = 5 is given in Sectional 

Theorem 1. For N = 0,1, 2, the function y( N \x) and the amplitude obtained by 
means of the above algorithm satisfy © up to the order 0(e N+l ) uniformly in x 6 [— 1, 1]: 

y (N)df2_ + ef{a (N) x)y (N) +x = 0{e N+ly (14) 
dx 

Moreover: 

(i) All of the coefficients y n ( x ) of y^ N \x) are continuous functions of x and satisfy y n (— 1) = 
y n (l) = and therefore y {N) (-l) = y {N) (l) = V N. 

(ii) If /3g( a o) 7^ (limit cycles of multiplicity higher than one are not considered here) then 
a2n+i = f2n+i{aQ,ai, ...,a2n+i',x) = for n = 0,1,2,.... The coefficients y 2 n+l(^) are odd 
functions of x and y2n( x ) are even functions of x. 
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(iii) When f{x) is a polynomial of degree q'mx 2 , then 

V2n{x) =p 2n (x 2 )(l-x 2 ) 3 / 2 , n = 1,2,3,..., 

?/2n+i(ic) = xp2n+i(a; 2 )(l - x 2 ), n = 0, 1,2, 
where p n (x) is a polynomial of degree nq — 1. 

Proof. If we replace y by y( N > and a by a^-* in @ and equate powers of e we obtain that 

y yo + x = (15) 

and, for n = 1,2, 3, every function y n (x) satisfies a first order linear differential equation 
which contains y , yi,...,y n _i: 

n n—1 

(a ,ai, ...,a m ;x)y n _ m _i — 0. (16) 

m=0 m=0 



The solution of (|T5|) which satisfies yo(— 1) = is obviously yo(x) = yl — x 2 . We solve 
recursively every one of the equations (|16j) for y n ( x ), n = 1,2, ...,iV, in terms of the 
preceding functions yo, yi,—, y n ~i and of the coefficients clq, a\,..., a n _i by imposing that 
y n (— 1) = for n = 1,2,3, ...,N. Then we obtain (|12|h Therefore, y( N ' satisfies (jHJ) with 
a replaced by up to the order 0(e N+1 ) and (|T1)) holds. 

We already have that y n (— 1) = for n = 0, 1, 2, iV. To show thesis (i), it remains 
to show that y n (l) = for n = 0, 1, 2, N when Eq. (|T3|) is assumed. We write (fT2]> in 
the form 

y n {x) = Wn (x) - /3 - l(a °' a . 1 7' a "- l) , (17) 

yo{x) 

with 

1 f 1 n— 1 n— 1 m /-l ~) 

w n (x) = pr < - J! 2/m(a:)y„_ m (x) + X X W / t k yn-m-i{t)f {k) {a t)dt \ (18) 

L / m=l m=0fc=0 • /s J 



and p n (ao,a%, ...,a n ) given in (|10|). Obviously, yo(x) = 0(y/l — x) when x — ► 1 and is a 
continuous function in [—1, 1]. From the hypothesis (|13j) . we have that /3 n (&o> o n ) = 

for n = 0, 1, 2, iV — 1 and then y n {x) = w n {x) for n = 1, 2, 3, N. From (|18|) . it is 
straightforward to show by induction over n that all of the functions w n {x) are continuous 
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in [—1, 1] and w n (x) = 0{y/l — x) when x — ► 1. Therefore, y n {x) is continuous in [—1, 1] 
and y n (l) = for n = 0, 1, 2, N. Then thesis (i) holds. 

From Definition 2, all of the odd functions /2n+i( a (b a l> •••> fl 2n+i; ^) vanish if the odd 
coefficients a 2n +\ vanish. This is so because it is necessary to have a couple (k, cr(k)) with 
odd k and <r(k) in the definition of b n ^ m in order to satisfy the condition Ylk=l ka(k) = 
2n + 1. Then, to show the first part of (ii) it suffices to show that a 2n +\ = 0, with 
n = 0,1,2,.... 

The function yo(x) is even and from (|12jl with n = 1: 

V\{x) = i-T f X f(a t)y {t)dt. (19) 

yo{x) J-i 

Taking into account that /3o( a o) = j\ f( a ot)yo{t)dt = and that f(x) is even, we 
have f (aot)yo(t)dt = ^/?o( a o) = 0. We can write the integral in ()19() in the form 
I-if(aot)y (t)dt = /° x f(a t)y (t)dt + Jq f(a t)y (t)dt = Jq f(a t)y (t)dt, which shows 
that y\{x) is odd. From (J10[) with n = 1 we have 

f3i(a , a{) = J ^ f (a x)y 1 (x)dx + a x J ^ xf'(a x)y (x)dx. 

The first integral vanishes because its integrand is odd and the second one equals Pb(ao), 
which is not zero from assumption. Therefore, /3i(ao, a\) = =^ a\ = and /i(ao, a>i;x) = 
0. Then, thesis (ii) is true for n = 0. Suppose that it is true for n = 0, 1,2, ...,m — 1 
with m > 0; we will show that thesis (ii) is also true for n = m. If thesis (ii) is true for 
n = 0, 1, 2, m — 1, then for n < m, the integrals in the first line of (|12j) may be written 
in the form 

L(n-l)/2j 



fc=0 
L("-1)/2J ,,- 



/x 
y n -2k-i(t)f2k(ao,ai, ■.■,a 2 k]t)dt = 

-!)/2J ,0 

Y] / J/n-2fc-i(*)/2fc(ao,ai,». ,a 2 k 4 ,t)dt + 

-)/2J 

Y] / yn-2k-i(t)f 2 k(ao,a 1 ,...,a 2 k;t)dt. (20) 

fc=0 70 

If f(x) is an even function of x, then all of the functions f 2 k(ao,ai, ...,a 2 k',x) are even. 
When n is even then y n -2k-l{P)i k = 0, 1, 2, [(» — 1) /2j , is odd and the function 
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f-iyn-2k-l{t)f2k( a 0, a l,---) a 2k'>t)dt are even functions of x. On the other hand, when n 
is odd then y n -2k-i(t), k = 0, 1, 2, \_(n — l)/2j , are even and the first sum in the right 
hand side of (|2Uj) equals /3 n _i(ao, a\, a n _i)/2 = and then the left hand side is an odd 
function of x. Obviously, the sum Y^kZi UkUn-k in <fl2|) is an even function of x if n < m 
is even and an odd function of x if n < m is odd. Then, the second part of thesis (ii) is 
true for n = m. 

From the induction hypothesis, a2 n +i = /2n+i = for n = 0,1,2,..., |/n/2j ~~ 1 an d 
then, for odd n = m, equation (jlOj) reads: 

)= / /m(ao, oi, a m ; x)yo(x)dx. (21) 

From Definition 2, for odd m, the function f m (ao, ai, a m ; x) is of the form / m (ao, ai, a m ; x) = 
x/ / (a x)a m +5 m (ai,a 3 , ...,a m _ 2 ;x), where g m (ai, a 3 , ...,a m _ 2 ;x) = when a x = a 3 = ... = 
o m _2 = 0. Then, for odd n = m equation ()21j) may be written in the form 

/?m(«0; Qlj •••) Om) = «m/3o( a o)- 

Therefore, the equation /3 m (ao, ai, a m ) = ^ a m = (and then / m = 0) for odd m. 
Then the first part of thesis (ii) is true for n = m. 

To show (hi) we recall the function y\ given in (|19|) . If f{x) is a polynomial of degree q 
in x 2 , then the integrand involved in the computation of y\{x) is a product of a polynomial 
of degree q in £ 2 by Vl -t 2 . An easy computation shows that 

/X 
f(aot)\/l — t 2 = a + xp(x 2 )\/\ — x 2 + (3 arcsin(x). 

where a and (5 are real numbers and p(x) a polynomial of degree q. But yi(±l) = and 
then 

J f{a t)Vl -t 2 = a ± /3arcsin(l) = 0^a = /3 = 0. 

Moreover, yi(x) is odd and from (|T§)) and Po( a o) = 0, yi(x) = 0(1 — x 2 ) when x — ► ±1. 
These facts mean that yi(x) = xpi(x 2 )(l — x 2 ), where p±(x) is a polynomial of degree 
q — 1. From these arguments and using induction over n in formula (|12j) we can easily 
obtain thesis (hi). 
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4 The first few terms of the expansion 



For e = 0, Eq. (j3J) has a continuum of periodic solutions: u(t) = ayo(t/a) for any value of 
a. Hence, yo(l) = does not impose any restriction over the parameter a. we determine 
do, o/v-i an d Ui(x),---,yN{x) by means of the above introduced algorithm. First, we take 



yo 0*0 = Vl — Then, for n = 0, 1,2, ...,iV — 1, solve the equation /3 n (ao, a%, o n ) = 
for a n and then compute We give below the details of this algorithm for N = 5. 



Definition 4. For convenience we define, for n = 2,4,6, the integrals 



7n(oo,a2, ...,a n _ 2 ) = J ^f(a x)y n (x)dx. 



The order re = . We calculate ao and 2/1(2) as follows. Compute the positive ao solutions 
of 



A)( a o) = / /(ooic)vi — £ 2 dx = 0. 



From (|12j) we obtain the expression (fT9|) for yi(x): 

1 



2/1 0*0 



yoO) 



f(a t)y (t)dt. 



The order re = 1 . We know that ai = and we calculate 2/2(2) from 1)121) with n = 2: 



2/2 (x) = - 



Vo(x) 



-ylix) + J* f(aot) yi (t)dt 



The order re = 2 . We calculate a 2 and 2/3(2). From (|10|) with re = 2 we have 

/? 2 (a ,0,a 2 ) = y ^ f(a x)y 2 (x)dx + a 2 J ^ xf'(a x)y (x)dx = 7 2 (a ) + a 2 /3 (a ) = 



Then, 

From l|12j) with re = 3: 
1 



a2 



72 (go) 
/5 («o)' 



(22) 



J/a(a:) 



2/0(2) 



2/1(2)2/2(2) + / ^f(a t)y 2 (t)dt + a 2 J tf(a o t)y (t)dt 
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The order n = 3 . We know that 03 = and we calculate 2/4 (x) from (|12|) with n = 4: 
1 



Va{x) 



yoO) 



yi{x)y 3 (x) + -yj(x) + J ^f{a t)y 3 (t)dt + a 2 J tf'(aot)yi(t)dt 



The order n = 4 . We calculate 04 and y${x). From (|10jl with n = 4 we have 

a 2 

/3 4 (ao,0,a 2 ,0,a 4 ) = 7 4 (a ,a 2 ) + a 2 7 2 (a ) + y/?o( a o) + a4/3 ( a o) = 



Then, 

From (|12j) with n = 5: 

, « 1 



74 (a , a 2 ) + a 2 7 2 (ao) + if#> (ao) , oq s 

° 4 = to) • (23) 



2/5 w 



yi(a;)y4(a:) + y2(x)y 3 (x) + / / ' (a Q t)y A {t)dt + a 2 / tf'(a t)y 2 (t)dt+ 



(24) 



i tf'(a t)y (t)dt + y J 1 t 2 f" (a t)y {t)dt 



Remark 1. The integral /?o( a o) coincides with the first Melnikov function or, equivalently, 
with the Abelian integral ^2] defined for the perturbed Hamiltonian system © whose 
level curves at e = are given by the set of circles x 2 + y 2 = a 2 . Thus, the equation 
/3o( a o) = is a nonlinear equation for ao with none, one or several solutions, whereas the 
remaining equations /3 n (ao, ai, a n ) = with n = 1,2, .., N — 1 are linear equations in a n 
with a unique solution for every ao solution of /3o(°o) = 0- Then, the number of positive 
solutions ao of /?o(°o) = determine the maximal number of limit cycles emerging from the 
period annulus when it is perturbed and the 0(1) approximation to their amplitude. The 
remaining equations /3 n ( a o> a l> ■••> a n) = with n = 1, 2, .., N — 1 determine recursively the 
coefficients ai,...,ajv-i, which are the pertubative correction to the first order amplitudes 
ao- In this way, for every solution of Po( a o) = 0, we completely determine the amplitude 

a and the form y{x) of the limit cycle solutions of (jHJ) up to the order 0(e ): 

N-l N 
a ~ a (N) = £ y{x) „ y (N) {x) = £ e n yn(x)> 

n=0 n=0 

If /(#) is a polynomial of degree g in x 2 , then /3o(^q) is a polinomial in a§ of degree q. 
Then, for small e, the maximum possible number of limit cycles of © is g, as it has been 
conjectured by Lins, Melo and Pugh (LMP-conjecture) 
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Remark 2. The approximate solution y( N \x) with amplitude aS^J has the same sym- 
metries than the exact solutions of Eq. ©: the inversion symmetry (x,y) <-> (— x, — y) 
and the parameter inversion symmetry (e,x,y) *-* (— e, —x,y). These symmetries are not 
destroyed at the perturbative level. 

Remark 3. From thesis (iii) of Theorem 1, when f(x) is a polynomial of degree q in x 2 
we have that 

LA72J W2J 



yW(x) = VT^ + x(l - x 2 ) £ P2n + i{x 2 )e 2n+1 + (1 - x 2 f 2 £ p 2n (x 2 )e 2n , 

n=0 n=l 

where p n (x) is a polynomial of degree nq — 1. 



5 Examples 

We perform the calculations proposed in Section 4 for two concrete examples, which are 
particular cases of the families 1 and 3 worked out in |13j . 

Example 1. The van der Pol oscillator is given for f(x) = x 2 — l. This system has a unique 
limit cycle, which is stable for e > 0. Hence, the only root of (3o(a) is a® = 2. For this 
value of ao, an approximate form of the limit cycle for small e is y^\x) = Y^ n =Q en Vn{x) 
with 



and 



2/0 0*0 = Vl - x 2 , yi(x) = x(l - x 2 ), 

y 2 (x) = ^(1 + 2x 2 )(l - x 2 ) 3 / 2 , y 3 (x) = ^x 3 (6x 2 - 5)(1 - x 2 ), 
y A (x) = ——(-4 - 23x 2 + 213x 4 - 156x 6 )(l - x 2 ) 3/2 



ys ( x ) = x 3 (10385 - 43794x 2 + 41424x 4 - 10560x 6 )(l - x 2 } 

yay ' 1036800 v 



An approximate form of its amplitude for small e is a(e) = + 0(e 8 ), with 

„ W = 2 + l e 2_Jg_ e 4 + 10M689 A 

96 552960 55738368000 v ' 
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We integrate Eq. © by a Runge-Kutta method in order to obtain the limit cycle. This 
curve is plotted in a continuous trace in Fig. l(a-b) for e = 0.8 and e = 2, respectively. 
The approximated limit cycle u^ 5 \x) = a^y^ (x/a^) is also plotted in those figures with 
a discontinuous trace for the same values of e. Let us remark that, in this case, even up 
to e = 3, the approximation u^(x) to the limit cycle is very good. Our calculation (j25j) 
agrees with the computational calculation of the 'exact' amplitudes given in Table 1, and 
also with the calculations on this system presented in [TO] . 

Example 2. The same process is performed for f(x) = 5x 4 — 9x 2 + 1. In this case, the 
system has two limit cycles, one stable and the other unstable. The polynomial Po(a) 



has two positive roots: ao = 




= 0.720677 (unstable limit cycle for e > 0) and 



a = = 1.75517 (stable limit cycle for e > 0). 

For the first limit cycle we have y^(x) = yo(%) + «/i(x) + e 2 y2(x) + e 3 ys(x) with 

y (x) = VI - x 2 , 

yi (x) = x(x 2 - 1) I 1 + x 2 j , 

y 2 ( x ) = (1 - x 2 ) 3/2 (0.06586 + 0.13173x 2 - 0.05241x 4 + 0.00505x 6 ) 

and 

y 3 (x) = x 3 (l - x 2 )(0.06081 - 0.11246x 2 + 0.05384x 4 - 0.00999x 6 + 0.00006x 8 ). 

An approximate expression of its amplitude for small e is 

o(e) = 0.720677 + 0.00390888 e 2 + C(e 4 ). 

For the second limit cycle we have y^ 3 \x) = yo(x) + ey\(x) + e 2 y2(x) + e 3 ys(x) with 

2/o (x) = VT - x 2 , 

f 2 ^(i 9^/41 + 61 2 \ 
y 1 {x) = x(x -1)11 — x I, 

y 2 (x) = (1 - x 2 ) 3/2 (0.98791 + 1.97583x 2 - 2.71374x 4 + 6.2545x 6 ) 
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and 



y 3 (x) = x 3 (l - x 2 )(0.846907 - 2.30045x 2 - 1.11829x 4 - 25.2371a; 6 + 28.2651x 8 ). 



An approximate form of its amplitude for small e is 

o(e) = 1.75517 - 0.0178803 e 2 + C(e 4 ). 



The comparison between the 'exact' and the approximated limit cycles can be seen in 
the plots of Fig. 2(a-b) for e = 0.8 and e = 2 and Fig. 3(a-b) for e = 0.4 and e = 0.8. 



6 Conclusions 



Periodic self-oscillations can arise in nonlinear systems. These are represented by isolated 
closed curves in phase space that we call limit cycles. The knowledge of the number, 
amplitude and shape of these solutions in a general nonlinear equation is an unsolved 
problem. 

In this work, a recursive algorithm to approximate the form y(x) and the amplitude 
a of limit cycles in the weakly nonlinear regime (e — > 0) of Lienard equations, x = y, 
y+ef(x)y + x = 0, with f(x) an even function, has been presented. The symmetries of the 
exact limit cycles are maintained at the perturbative level. Several examples showing the 
application of this scheme have been given, and the accuracy of the method in comparison 
with the direct numerical integration has also been tested. In particular, we have explicitly 
obtained the C(e 8 ) approximation for the amplitude of the van der Pol limit cycle. 
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Table 1. The value ax represents the approximated amplitude a(e) of the van der Pol 
limit cycle obtained from (|25|) for the indicated values of e. The value aE represents the 
amplitude a obtained by integrating directly the system with a Runge-Kutta method. 



e 


0.1 


0.2 


0.3 


0.4 


0.5 


ar 


2.000104 


2.000414 


2.000922 


2.001619 


2.002488 


a E 


2.00010 


2.00041 


2.00092 


2.00161 


2.00248 


€ 


0.6 


0.7 


0.8 


0.9 


1 




2.003509 


2.004658 


2.005906 


2.007222 


2.008569 


a E 


2.00351 


2.00466 


2.00591 


2.00724 


2.00862 



Figures 

Fig. la-b: Exact limit cycle (continuous line) and approximated limit cycle (dis- 
continuous line) up to order 0(e 5 ) for the van der Pol system with (a) e = 0.8 and (b) 
e = 2. 

Fig. 2a-b: Exact limit cycle (continuous line) and approximated limit cycle (dis- 
continuous line) up to order C(e 3 ) for the unstable limit cycle of the example 2 with (a) 
e = -0.8 and (b) e = -2. 

Fig. 3a-b: Exact limit cycle (continuous line) and approximated limit cycle (discon- 
tinuous line) up to order 0(e 3 ) for the stable limit cycle of the example 2 with (a) e = 0.4 
and (b) e = 0.8. 
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